function [min_scale,max_scale,min_scalem,max_scalem,min_scalef,max_scalef] = get_scale(singlefemale,singlemale,couple,leisurem,leisuref, hworkm, hworkf, consumption,wagem,wagef,nonlabor,divorce_costsM,divorce_costsF,divorce_costsMF,ub)
% Program uses YALMIP (https://yalmip.github.io) connected to CPLEX or GUROBI

% INPUT VARIABLES:
% leisurem[households]      -- leisure per week of the male
% leisuref[households]      -- leisure per week of the female
% consumption[households]   -- weekly non-assignable consumption
% wagem[households]         -- wage rate of the male
% wagef[households]         -- wage rate of the female
% hworkm[households]        -- time spend on home production by the male
% hworkf[households]        -- time spend on home production by the female
% nonlabor[households]      -- weekly nonlabor income of the household


%defining constants:
workh = 112;        % total potential working hours
num_of_couples = length(find(couple)); %number of couples in the data set
num_of_hh = length(couple); %total number of households in the data.
e = 0.0000001;

%% DEFINING DECISION VARIABLES
privateqM = sdpvar(num_of_hh,1);
privateqF = sdpvar(num_of_hh,1);
privatemhworkm = sdpvar(num_of_hh,1);
privatemhworkf = sdpvar(num_of_hh,1);
privatefhworkm = sdpvar(num_of_hh,1);
privatefhworkf = sdpvar(num_of_hh,1);
privateleisurem = sdpvar(num_of_hh,1);
privateleisuref = sdpvar(num_of_hh,1);

scale = sdpvar(num_of_couples,1);
scalem = sdpvar(num_of_couples,1);
scalef = sdpvar(num_of_couples,1);

aleisurem = sdpvar(1);
aleisuref = sdpvar(1);
ahworkm = sdpvar(1);
ahworkf = sdpvar(1);
aq = sdpvar(1);

nl_incomeM = sdpvar(num_of_hh,1);
nl_incomeF = sdpvar(num_of_hh,1);

%% SETTING UP THE CONSTRAINTS
Constraints = [];

% private consumption should be balanced %

% Hicksian good
for i=1:num_of_hh
    Constraints=[Constraints, ...
        consumption(i)*(1-aq) == privateqM(i) + privateqF(i), privateqM(i) >=0, privateqF(i) >= 0];
end
for i=1:num_of_hh
      if singlefemale(i) == 1
            Constraints = [Constraints, ...
                privateqM(i) == 0];
        else
            if singlemale(i) == 1
               Constraints = [Constraints,...
                  privateqF(i) == 0];
            end
      end 
end

% Domestic work by male 
for i=1:num_of_hh
    Constraints=[Constraints, ...
        hworkm(i)*(1-ahworkm) == privatemhworkm(i) + privatemhworkf(i), privatemhworkm(i) >= 0, privatemhworkf(i) >= 0];
end
for i=1:num_of_hh
      if singlefemale(i) == 1
            Constraints = [Constraints, ...
                privatemhworkm(i) == 0];
        else
            if singlemale(i) == 1
               Constraints = [Constraints,...
                  privatemhworkf(i) == 0];
            end
      end 
end

% Domestic work by female
for i=1:num_of_hh
    Constraints=[Constraints, ...
        hworkf(i)*(1-ahworkf) == privatefhworkm(i) + privatefhworkf(i), privatefhworkm(i) >=0 , privatefhworkf(i) >= 0];
end
for i=1:num_of_hh
      if singlefemale(i) == 1
            Constraints = [Constraints, ...
                privatefhworkm(i) == 0];
        else
            if singlemale(i) == 1
               Constraints = [Constraints,...
                  privatefhworkf(i) == 0];
            end
      end 
end

% Leisure male
for i=1:num_of_hh
    Constraints=[Constraints, ...
        leisurem(i)*(1-aleisurem) == privateleisurem(i)];
end

% Leisure female
for i=1:num_of_hh
    Constraints=[Constraints, ...
        leisuref(i)*(1-aleisuref) == privateleisuref(i)];
end

% balanced nonlabor income:
for i=1:num_of_hh
        if couple(i) == 1
            Constraints = [Constraints, nonlabor(i) == nl_incomeM(i)+nl_incomeF(i)];
        else
            if singlefemale(i) == 1
                Constraints = [Constraints, nonlabor(i) == nl_incomeF(i)];
            else
                if singlemale(i) == 1
                Constraints = [Constraints, nonlabor(i) == nl_incomeM(i)];
                end
            end 
        end
end

% bounding the possible sharing of nonlabor income between partners [.4,.6]
for i=1:num_of_hh
    % check that dealing with couple
    if couple(i) == 1
        if nonlabor(i)>=0
            Constraints = [Constraints,
                .4*nonlabor(i) <= nl_incomeM(i) <= .6*nonlabor(i), .4*nonlabor(i) <= nl_incomeF(i) <= .6*nonlabor(i)];
        else 
            Constraints = [Constraints,
                .6*nonlabor(i) <= nl_incomeM(i) <= .4*nonlabor(i), .6*nonlabor(i) <= nl_incomeF(i) <= .4*nonlabor(i)];
        end
    end
end

% INDIVIDUAL RATIONALITY:
for i=1:num_of_couples
        % IR-M
        Constraints = [Constraints,
            (wagem(i)*workh)*(divorce_costsM(i)-e) + nl_incomeM(i) <= privateleisurem(i)*wagem(i) + privatemhworkm(i)*wagem(i) + privatefhworkm(i)*wagef(i) + privateqM(i) + hworkm(i)*ahworkm*wagem(i) + hworkf(i)*ahworkf*wagef(i) + aq*consumption(i) + leisurem(i)*aleisurem*wagem(i) + leisuref(i)*aleisuref*wagef(i)];
        % IR-W
        Constraints = [Constraints,
            (wagef(i)*workh)*(divorce_costsF(i)-e) + nl_incomeF(i) <= privateleisuref(i)*wagef(i) + privatemhworkf(i)*wagem(i) + privatefhworkf(i)*wagef(i) + privateqF(i) + hworkm(i)*ahworkm*wagem(i) + hworkf(i)*ahworkf*wagef(i) + aq*consumption(i) + leisurem(i)*aleisurem*wagem(i) + leisuref(i)*aleisuref*wagef(i)];
end

% NO BLOCKING PAIR CONDITIONS
for i=1:num_of_couples
    for j=1:num_of_hh
        if (couple(i)== 1) && (couple(j)==1)
            Constraints = [Constraints, 
                (wagem(i) + wagef(j))*workh*(divorce_costsMF(i,j)-e) + nl_incomeM(i) + nl_incomeF(j) <= privateleisurem(i)*wagem(i) + privatemhworkm(i)*wagem(i) + privatefhworkm(i)*wagef(j) + privateqM(i) + privateleisuref(j)*wagef(j) + privatemhworkf(j)*wagem(i) + privatefhworkf(j)*wagef(j) + privateqF(j) + max(hworkm(i),hworkm(j))*ahworkm*wagem(i) + max(hworkf(i),hworkf(j))*ahworkf*wagef(j) + aq*max(consumption(i),consumption(j)) + max(leisurem(i),leisurem(j))*aleisurem*wagem(i) + max(leisuref(i),leisuref(j))*aleisuref*wagef(j)];
        end
        if (couple(i)==1) && (singlefemale(j)==1)
           Constraints = [Constraints, 
           (wagem(i) + wagef(j))*workh*(divorce_costsMF(i,j)-e) + nl_incomeM(i) + nl_incomeF(j) <= privateleisurem(i)*wagem(i) + privatemhworkm(i)*wagem(i) + privatefhworkm(i)*wagef(j) + privateqM(i) + privateleisuref(j)*wagef(j) + privatemhworkf(j)*wagem(i) + privatefhworkf(j)*wagef(j) + privateqF(j) + max(hworkm(i),hworkm(j))*ahworkm*wagem(i) + max(hworkf(i),hworkf(j))*ahworkf*wagef(j) + aq*max(consumption(i),consumption(j)) + max(leisurem(i),leisurem(j))*aleisurem*wagem(i) + max(leisuref(i),leisuref(j))*aleisuref*wagef(j)];
        end
       if (couple(i)==1) && (singlemale(j)==1)
                    Constraints = [Constraints, 
           (wagem(j) + wagef(i))*workh*(divorce_costsMF(j,i)-e) + nl_incomeM(j) + nl_incomeF(i) <= privateleisurem(j)*wagem(j) + privatemhworkm(j)*wagem(j) + privatefhworkm(j)*wagef(i) + privateqM(j) + privateleisuref(i)*wagef(i) + privatemhworkf(i)*wagem(j) + privatefhworkf(i)*wagef(i) + privateqF(i) + max(hworkm(j),hworkm(i))*ahworkm*wagem(j) + max(hworkf(j),hworkf(i))*ahworkf*wagef(i) + aq*max(consumption(j),consumption(i)) + max(leisurem(j),leisurem(i))*aleisurem*wagem(j) + max(leisuref(j),leisuref(i))*aleisuref*wagef(i)];
        end
    end
end
   
Constraints = [Constraints,  0 <= aq <= 1, 0 <= aleisurem <= ub, 0 <= aleisuref <= ub, 0 <= ahworkm <= 1, 0 <= ahworkf <= 1];

for i=1:num_of_hh
    if (couple(i)==1)
    Constraints = [Constraints,...
        scale(i) == ((1+aleisurem)*wagem(i)*leisurem(i) + (1+aleisuref)*wagef(i)*leisuref(i) + (1+ahworkm)*wagem(i)*hworkm(i)+ (1+ahworkf)*wagef(i)*hworkf(i) + (1+aq)*consumption(i))/(wagem(i)*workh + wagef(i)*workh + nonlabor(i)),...
        scalem(i) == (privateleisurem(i)*wagem(i) + privatemhworkm(i)*wagem(i) + privatefhworkm(i)*wagef(i) + privateqM(i) + aleisurem*wagem(i)*leisurem(i) +  aleisuref*wagef(i)*leisuref(i) + ahworkm*wagem(i)*hworkm(i)+ ahworkf*wagef(i)*hworkf(i) + aq*consumption(i))/(wagem(i)*workh + wagef(i)*workh + nonlabor(i)),...
        scalef(i) == (privateleisuref(i)*wagef(i) + privatemhworkf(i)*wagem(i) + privatefhworkf(i)*wagef(i) + privateqF(i) + aleisurem*wagem(i)*leisurem(i) +  aleisuref*wagef(i)*leisuref(i) + ahworkm*wagem(i)*hworkm(i)+ ahworkf*wagef(i)*hworkf(i) + aq*consumption(i))/(wagem(i)*workh + wagef(i)*workh + nonlabor(i))];
    end
end


%% Specifying Objective Function
u = sdpvar(num_of_couples,1);    % vector of parameters necessasry to specify the function
Objective1 = sum(u.*scale);          % objective function that uses parameter -- u would be using as +/- unit vector 
Objective2 = sum(u.*scalem);          % objective function that uses parameter -- u would be using as +/- unit vector 
Objective3 = sum(u.*scalef);          % objective function that uses parameter -- u would be using as +/- unit vector 


% Solving the problem
options = sdpsettings('verbose',0,'solver','gurobi');
Opt1 = optimizer(Constraints,Objective1,options,u,scale); 
Opt2 = optimizer(Constraints,Objective2,options,u,scalem); 
Opt3 = optimizer(Constraints,Objective3,options,u,scalef); 


% parametrizations of the problem
param = eye(num_of_couples,num_of_couples);

for i=1:num_of_couples
    % economies of scale
    [xvalue, errorcode] = Opt1(param(:,i));
    if errorcode~=0
        display('Something went wrong (min)');    
    end
    min_scale(i) = xvalue(i);
       
    [xvalue, errorcode] = Opt1(-param(:,i));
    if errorcode~=0
        display('Something went wrong (max)');
    end
    max_scale(i) = xvalue(i);
    
    % male RICEB
    [xvalue, errorcode] = Opt2(param(:,i));
    if errorcode~=0
        display('Something went wrong (min)');
    end
    
    min_scalem(i) = xvalue(i);
    [xvalue, errorcode] = Opt2(-param(:,i));
    if errorcode~=0
        display('Something went wrong (max)');
    end
    max_scalem(i) = xvalue(i);
    
    % female RICEB
    [xvalue, errorcode] = Opt3(param(:,i));
    if errorcode~=0
        display('Something went wrong (min)');
    end
    
    min_scalef(i) = xvalue(i);
    [xvalue, errorcode] = Opt3(-param(:,i));
    if errorcode~=0
        display('Something went wrong (max)');
    end
    max_scalef(i) = xvalue(i);
   
end



return